The study domain and the model domain setup
#
fontsize_tmp = fontsize
fontsize=20
from mpl_toolkits.axes_grid1.inset_locator import zoomed_inset_axes, inset_axes
from mpl_toolkits.axes_grid1.inset_locator import mark_inset
topocmap = plt.cm.gist_earth
topocmap.set_bad('w', alpha=0)
#First the greater area map
with nc(os.path.join(geodataf,'topo_2.nc')) as topof:
Gtopo=np.ma.masked_equal(topof.variables['topo'][:],0)
Gtopolon=topof.variables['longitude'][:]
Gtopolat=topof.variables['latitude'][:]
#Now only the tiwi-islands manp
with nc(os.path.join(geodataf,'Topo.nc')) as topof:
topo=np.ma.masked_equal(topof.variables['topo'][:],0)
topolon=topof.variables['longitude'][:]
topolat=topof.variables['latitude'][:]
with nc(CPOLF) as fnc:
lon=fnc.variables['longitude'][:]
lat=fnc.variables['latitude'][:]
fig = plt.figure(figsize=(20,20))
ax = fig.add_subplot(111)
ls = LightSource(azdeg=315, altdeg=45)
M = Basemap(projection='rotpole', llcrnrlat=min(Gtopolat), llcrnrlon=min(Gtopolon), urcrnrlat=max(Gtopolat), urcrnrlon=max(Gtopolon),
resolution='c', area_thresh=1, ax=ax, lon_0=0, o_lon_p=0., o_lat_p=90.,)
borders = (([122.963, 139.055], [-17.091, -7.443]),
([127.385, 134.693], [-14.505, -10.005]),
([129.557, 132.529], [-13.757, -10.757]))
text=(('UM 4km', fontsize+4), ('UM 1.33km', fontsize+4), ('UM 0.44km', fontsize+4))
M.pcolormesh(Gtopolon, Gtopolat, ls.hillshade(Gtopo[:], vert_exag=1), cmap= topocmap, ax=ax)
#M.drawmapscale(M.boundarylons[0]-0.7, M.boundarylats[0]+0.7, max(topolon), min(topolat), 10,
# barstyle='fancy', fontsize=19, labelstyle='simple')
for conf in range(3):
lons,lats = borders[conf]
txt, fts = text[conf]
x, y = M([lons[0], lons[0], lons[1],lons[1]], [lats[0], lats[1], lats[1], lats[0]])
xy = list(zip(x,y))
ttt = plt.text(*M(lons[0]+0.1, lats[0]+0.05), txt, fontsize=fts)
ttt.set_bbox(dict(facecolor='w', alpha=0.7, edgecolor='none', lw=1))
poly = Polygon(xy, alpha=1, edgecolor='k', lw=1, facecolor='none')
plt.gca().add_patch(poly)
#M.drawcoastlines()
M.scatter(*M(131.043101, -12.250323), marker='*',color='k', s=200)
plt.text(*M(borders[0][0][0]+0.15, borders[0][-1][1]-0.6), 'b)', fontsize=fontsize+12)
axins = zoomed_inset_axes(ax, 7, loc=2)
axins.set_xlim(min(topolon), max(topolon))
axins.set_ylim(min(topolat), max(topolat))
tmap = Basemap(projection='cyl',llcrnrlat=min(lat), llcrnrlon=min(lon), urcrnrlat=max(lat), urcrnrlon=max(lon),
resolution='c', area_thresh=1, ax=axins)
tmap.pcolormesh(topolon, topolat, ls.hillshade(topo[:], vert_exag=1), cmap= topocmap)
mark_inset(ax, axins, loc1=3, loc2=4, fc="k", ec='k', lw=0.5)
plt.annotate('a)', xy=(0.01 ,0.92), xycoords='axes fraction', fontsize=fontsize+12)
lons, lats = [min(Gtopolon), max(Gtopolon)], [min(Gtopolat), max(Gtopolat)]
axins2 = inset_axes(ax, width='20%', height='20%', loc=4)
Worldmap= Basemap(projection='ortho', lat_0=Gtopolat[int(len(Gtopolat)/2)], lon_0=Gtopolon[int(len(Gtopolon)/2)],
ax=axins2)
Worldmap.drawcoastlines(color='k', linewidth=0.7)
Worldmap.fillcontinents(color='gray')
bx, by = Worldmap(M.boundarylons, M.boundarylats)
xy = list(zip(bx,by))
poly = Polygon(xy, alpha=1, edgecolor='none', lw=0, facecolor='darkcyan')
_ = Worldmap.ax.add_patch(poly)
_ = fig.subplots_adjust(bottom=0.35)
_ = fig.savefig('Vid/Topo.png', bbox_inches='tight', format='png', dpi=300)
fontsize = fontsize_tmp
Now calculate the map of spectral variance in the diurnal cycle
Burst vs break
# Define extreme Events in the CPOL era
extr1h = pd.read_hdf(extremeTS, 'P10min')
monsoon10m = pd.read_pickle(BurstF)
monsoon10m = monsoon10m.dropna()
monsoon1h = monsoon10m.groupby(pd.TimeGrouper('1h')).mean()
breaks = monsoon10m.loc[monsoon10m == 0].index
bursts = monsoon10m.loc[monsoon10m > 0].index
thresh = perc['10min'][95]
S = pd.Series(np.arange(len(times)),index=times)
print('Bursts: %0.2f %% Breaks: %0.2f %%'
%(len(bursts)/len(monsoon10m.index) * 100, len(breaks)/len(monsoon10m.index) * 100 ))
Bursts: 39.73 % Breaks: 60.27 %
#Plot the diurnal Cycle
%matplotlib inline
nplot=1
import matplotlib.gridspec as gridspec
gs = gridspec.GridSpec(5, 3)
sns.set_style("whitegrid", {'axes.grid':True, 'ticks':True})
from matplotlib.ticker import ScalarFormatter
#Plot the output map
fig = plt.figure(figsize=(15,8))
cmap = colmap2.Blues
cmap2 = colmap.GMT_drywet
cmap.set_under('w'),cmap.set_bad('w')
cmap2.set_under('w'),cmap2.set_bad('w')
ax1 = plt.subplot(gs[:2, 0])
#ax1 = fig.add_subplot(231)
ax1.set_title('Rainfall Climatology', size=fontsize)
X, Y = np.meshgrid(lon.filled(-50), lat.filled(-50))
im = m1.contourf(X, Y, map1, np.arange(0,13,2), cmap=colmap.GMT_drywet, ax=ax1)
#im = m1.pcolormesh(lon.filled(-50), lat.filled(-50), map1, vmin=0, cmap=colmap.GMT_drywet, ax=ax1)
m1.drawcoastlines(linewidth=.7)
cbar = m1.colorbar(im , location='bottom', pad="5%")
cbar.set_label('Rain-Rate [mm/d]', size=fontsize)
cbar.ax.tick_params(labelsize=fontsize-4)
sns.set_style("darkgrid", {'axes.grid':True, 'ticks':True})
ax2 = plt.subplot(gs[:2, 1])
sns.lineplot('x', "y", data=sns_data, ax=ax2, lw=1)
ax2.set(xscale="log")
ax2.set_xlim(0.25,60)
ax2.set_ylim(0,0.6)
ax2.set_xlabel('Periods ($\\tau$) [d]',size=fontsize)
ax2.set_title('Spectral variance [\%]',size=fontsize)
ax2.tick_params(labelsize=fontsize-4)
ax2.set_ylabel('')
ax2.xaxis.set_major_formatter(ScalarFormatter())
sns.set_style("whitegrid", {'axes.grid':True, 'ticks':True})
ax3 = plt.subplot(gs[:2, 2])
im2 = m1.contourf(X, Y, map2, np.arange(0, 9)/10. ,cmap='Blues', ax=ax3)
m1.drawcoastlines(linewidth=.7)
ax3.set_title('Diurnal Cycle', size=fontsize)
cbar = m1.colorbar(im2,location='bottom',pad="5%")
cbar.set_label('Spectral variance [\%]', size=fontsize)
cbar.ax.tick_params(labelsize=fontsize-4)
# %%
#fig.subplots_adjust(right=0.98, bottom=0.03, top=0.99,left=0.01, hspace=0, wspace=0.17)
#'''
sns.set_style("darkgrid", {'axes.grid':True, 'ticks':True})
#fig, ax = plt.subplots(1, figsize=(15,10))
ax4 = plt.subplot(gs[3:, :])
sns.lineplot(x='Local Time', y='Rain-Rate', data=diurnacl_cycle, hue='Monsoon-Stage', lw=3, ax=ax4)
ax4.set_ylabel('Rain-Rate mm/d', fontsize=fontsize)
ax4.set_xlabel('Local Time [h]', fontsize=fontsize)
#ax.xaxis.set_major_formatter(matplotlib.dates.DateFormatter('%_H'))
ax4.legend(loc=0, fontsize=fontsize-4)
ax4.tick_params(labelsize=fontsize)
fig.subplots_adjust(bottom=0.01, top=0.98, left=0.01, right=0.98, hspace=0.1, wspace=0.2)
#fig.subplots_adjust(bottom=0.145, top=0.99, right=0.99, left=0.1)
nplot += 1; _ = fig.savefig('Vid/Fig_%03i.png'%nplot, bbox_inches='tight', format='png', dpi=300)
The above plots show that diurnal cycle is the most important mode of variability in the considered area. The island thunderstrom hector is the most extreme convective system that occurs at the time-scale at perdiods of 24 h. The Hector 'signal' is most pronounced during the Monsoon break period. Before studying Hector we split the CPOL area into break and burst period and investigate the occurrence of extremes during Monsoon break and burst.
First we define extreme events by applying the 99th percentile threshold of the total 10min dataset. Then we define the break and burst events by the definition of Narsey et al. and calculate the number of burst and the number of break events during the CPOL era
# Now Plot the maps
mpld3.disable_notebook()
sns.set_style("whitegrid", {'axes.grid':True, 'ticks':True})
maxval =[(0,15), (0,15), (-12,12), (0, 150), (0, 150), (-200, 200)]
names=['Breaks', 'Bursts', 'Breaks - Bursts',
'Extremes During Breaks', 'Extremes During Bursts', 'Breaks - Bursts']
cmap_list = [(colmap.GMT_drywet,'k'), (colmap.GMT_drywet,'k'), (colmap.GMT_polar_r,'k'),
(mpl.cm.magma,'w'), (mpl.cm.magma,'w'), (colmap.GMT_polar_r,'k')]
C = [Mean[1]*1.5, Mean[2], Mean[1]*1.5 - Mean[2], D[1]*1.2, D[2], D[1]*1.2 - D[2]]
fig = plt.figure(figsize=(15,10))
mask = ((D[0] * 0 )+1)
im=[]
for i, idx in enumerate(C):
ax = fig.add_subplot(2,3,i+1)
data = np.zeros([len(lat),len(lon)])
colm, coast = cmap_list[i]
colm.set_bad(dict(k='w', w='k')[coast])
colm.set_under(dict(k='w', w='k')[coast])
try:
im.append(m.pcolormesh(lon, lat, mask*C[i].filled(0),vmin=maxval[i][0],vmax=maxval[i][1],cmap=colm))
except AttributeError:
im.append(m.pcolormesh(lon, lat, mask*C[i],vmin=maxval[i][0],vmax=maxval[i][1], cmap=colm))
m.drawcoastlines(color=coast)
ax.set_title('%s'%names[i],size=fontsize)
fig.subplots_adjust(right=0.98, bottom=0.01, top=0.98,left=0.01, hspace=0, wspace=0)
cbar_ax1 = fig.add_axes([0.05, 0.57,0.58, 0.015])
cbar1 = fig.colorbar(im[0], cax=cbar_ax1, orientation='horizontal')
cbar1.ax.tick_params(labelsize=fontsize-4)
cbar1.set_label(u'Rain-Rate [mm/d]',size=fontsize)
cbar_ax2 = fig.add_axes([0.685, 0.57,0.27, 0.015])
cbar2 = fig.colorbar(im[2], cax=cbar_ax2, orientation='horizontal')
cbar2.ax.tick_params(labelsize=fontsize-4)
cbar2.set_label(u'Difference',size=fontsize)
cbar_ax3 = fig.add_axes([0.05, 0.09,0.58, 0.015])
cbar3 = fig.colorbar(im[3], cax=cbar_ax3, orientation='horizontal')
cbar3.ax.tick_params(labelsize=fontsize-4)
cbar3.set_label(u'N. of events above 99th%',size=fontsize)
cbar_ax4 = fig.add_axes([0.685, 0.09,0.27, 0.015])
cbar4 = fig.colorbar(im[-1], cax=cbar_ax4, orientation='horizontal')
cbar4.ax.tick_params(labelsize=fontsize-4)
cbar4.set_ticks([-200, -100, 0, 100, 200])
cbar4.set_label(u'Difference',size=fontsize)
nplot+=1 ; _ = fig.savefig('Vid/Fig_%03i.png'%nplot, bbox_inches='tight', format='png', dpi=300)
Now we compare the ensemble timeseries of the simulated rainfall for the 1.33km and 0.44km Simulations with the CPOL observations.
time, obs, rain = [], [], []
for e in range(len(ensembles)):
r1 = list(UM133ens.variables['lsrain'][e].mean(axis=(1,2,3)).values)
r2 = list(UM044ens.variables['lsrain'][e].mean(axis=(1,2,3)).values)
obs+=len(r1)*['UM 1.33km'] + len(r2)*['UM 0.44km']
t1 = pd.DatetimeIndex(UM133ens.coords['t'].values).tz_localize(utc).tz_convert(timezone).tz_localize(None).to_pydatetime()
rain += r1 + r2
t2 = pd.DatetimeIndex(UM044ens.coords['t'].values).tz_localize(utc).tz_convert(timezone).tz_localize(None).to_pydatetime()
time += list(t1) + list(t2)
r3 = list(OBS.dataset[1].variables['lsrain'].mean(axis=(1,2)).values)
t3 = list(pd.DatetimeIndex(OBS.dataset[1].coords['t'].values).tz_localize(utc).tz_convert(timezone).tz_localize(None).to_pydatetime())
o3 = len(r3)*['CPOL']
rain += r3
obs += o3
time += t3
plot_data = pd.DataFrame({'rain-rate':rain, 'time':time, 'Type':obs})
#Plot area avg TS
import matplotlib.dates as mdates
mpld3.disable_notebook()
sns.set_style("darkgrid", {'axes.grid':True, 'ticks':True})
myFmt = mdates.DateFormatter('%d/%m')
#
#Get the minimum time period that is covered by all simulations and also observations
fig = plt.figure(figsize=(15,7))
ax = fig.add_subplot(111)
sns.lineplot(x='time', y='rain-rate', hue='Type', ax=ax, data=plot_data, lw=1.3)
plt.subplots_adjust(left=0.1, bottom=0.1, right=0.98, top=0.9)
ax.legend(loc=0, fontsize=fontsize-4)
ax.tick_params(labelsize=fontsize-4)
#ax.set_xlim(plot_data.time.iloc[0], plot_data.time.iloc[-1])
ax.set_ylim(-0.001,0.5)
ax.set_ylabel('Rain-Rate [mm/h]', fontsize=fontsize)
_ = ax.xaxis.set_major_formatter(myFmt)
_ = ax.set_xlabel('')
#_ = ax.set_title('Comparison of Area-Avg Rainfall', fontsize=fontsize)
fig.subplots_adjust(bottom=0.145, top=0.9, right=0.99, left=0.1)
nplot+=1; _ = fig.savefig('Vid/Fig_%03i.png'%nplot, bbox_inches='tight', format='png', dpi=300)
Now lets look at some animations of UM rainfall occuring over the Tiwi-islands and compare it to the cpol observations
os.chdir('/home/unimelb.edu.au/mbergemann/Work/Programming/Extremes/python/')
#Plot Animatoin
import io
import base64
from IPython.display import HTML
outvid=os.path.join('Vid','WeekOfHector-Ens-3.mp4')
video = io.open(outvid, 'r+b').read()
encoded = base64.b64encode(video)
HTML(data='''<video alt="test" width="950" height="500" loop="true" controls>
<source src="data:video/mp4;base64,{0}" type="video/mp4" />
</video>'''.format(encoded.decode('ascii')))
#Create the avg. maps of rainfall 0.44km
#m1, m2 = None, None
mpld3.disable_notebook()
sns.set_style("whitegrid", {'axes.grid':True, 'ticks':True})
rom = {1:'i', 2:'ii', 3:'iii', 4:'iv', 5:'v', 6:'vi', 7:'vii', 8:'viii'}
fig = plt.figure(figsize=(15,12))
UM = UM044ens.groupby('t.day').sum(dim='t').mean(dim='day')
obs = OBS.dataset[1].groupby('t.day').sum(dim='t').mean(dim='day')
colm = colmap2.Blues
colm.set_under('w', alpha=0)
colm.set_bad('w', alpha=0)
tsteps = pd.DatetimeIndex(OBS.dataset[1].coords['t'].values).tz_localize(utc).tz_convert(timezone).to_pydatetime()
obs_data = np.ma.masked_less(obs.variables['lsrain'][:].values, -0.01) / 6
um_data = np.ma.masked_less(UM.variables['lsrain'][:].values, -0.01) / 6
lat1, lon1 = obs.variables['latitude'][:,0], obs.variables['longitude'][0,:]
lat2, lon2 = UM.coords['lat'][:], UM.coords['lon'][:]
ax = [fig.add_subplot(3,3,1)]
if m2 is None:
m2 = Basemap(llcrnrlat=min(lat2), llcrnrlon=min(lon2), urcrnrlat=max(lat2), urcrnrlon=max(lon2),
resolution='f', area_thresh=1)
if m1 is None:
m1 = Basemap(llcrnrlat=min(lat1), llcrnrlon=min(lon1), urcrnrlat=max(lat1), urcrnrlon=max(lon1),
resolution='f', area_thresh=1)
#m[0].pcolormesh(topolon, topolat, ls.hillshade(topo[:], vert_exag=1), cmap='gray', ax=ax[0])
im = [m2.pcolormesh(lon1, lat1, obs_data,vmin=0.0,vmax=5,cmap=colm)]
m2.drawcoastlines(linewidth=0.8)
cbar_ax = fig.add_axes([0.12, 0.12, 0.74, 0.02])
cbar=fig.colorbar(im[-1], cax=cbar_ax, orientation='horizontal')
cbar.ax.tick_params(labelsize=24)
cbar.set_label('Rain-Rate [mm/d]',size=24)
tit = ['CPOL']
ax[-1].set_title('CPOL', fontsize=fontsize)
for i in range(1,9):
if m2 is None:
m2 = Basemap(llcrnrlat=min(lat2), llcrnrlon=min(lon2), urcrnrlat=max(lat2), urcrnrlon=max(lon2),
resolution='f', area_thresh=1)
ax.append(fig.add_subplot(3,3,i+1))
m2.drawcoastlines(linewidth=0.8)
im.append(m2.pcolormesh(lon2, lat2, um_data[i-1][0],vmin=0.0,vmax=5,cmap=colm))
tit.append('Ens. %s'%rom[i].upper())
ax[-1].set_title(tit[-1],fontsize=fontsize)
fig.subplots_adjust(right=0.99, bottom=0.15, top=0.95,left=0.01, hspace=0.2, wspace=0.0)
nplot+=1 ; _ = fig.savefig('Vid/Fig_%03i.png'%nplot, bbox_inches='tight', format='png', dpi=300)
#Create the avg maps of rainfall for 1.33km
sns.set_style("whitegrid", {'axes.grid':True, 'ticks':True})
fig = plt.figure(figsize=(15,12))
UM = UM133ens.groupby('t.day').sum(dim='t').mean(dim='day')
obs = OBS.dataset[1].groupby('t.day').sum(dim='t').mean(dim='day')
colm = colmap2.Blues
colm.set_under('w', alpha=0)
colm.set_bad('w', alpha=0)
tsteps = pd.DatetimeIndex(OBS.dataset[1].coords['t'].values).tz_localize(utc).tz_convert(timezone).to_pydatetime()
obs_data = np.ma.masked_less(obs.variables['lsrain'][:].values, -0.01) / 6
um_data = np.ma.masked_less(UM.variables['lsrain'][:].values, -0.01) / 6
lat1, lon1 = obs.variables['latitude'][:,0], obs.variables['longitude'][0,:]
lat2, lon2 = UM.coords['lat'][:], UM.coords['lon'][:]
ax = [fig.add_subplot(3,3,1)]
if m1 is None:
m1 = Basemap(llcrnrlat=min(lat1), llcrnrlon=min(lon1), urcrnrlat=max(lat1), urcrnrlon=max(lon1), resolution='f',
area_thresh=1)
#m[0].pcolormesh(topolon, topolat, ls.hillshade(topo[:], vert_exag=1), cmap='gray', ax=ax[0])
im = [m2.pcolormesh(lon1, lat1, obs_data,vmin=0.0,vmax=5,cmap=colm)]
m2.drawcoastlines(linewidth=0.8)
cbar_ax = fig.add_axes([0.12, 0.12, 0.74, 0.02])
cbar=fig.colorbar(im[-1], cax=cbar_ax, orientation='horizontal')
cbar.ax.tick_params(labelsize=24)
cbar.set_label('Avg. Rain-Rate [mm/d]',size=24)
tit = ['CPOL']
ax[-1].set_title('CPOL', fontsize=fontsize)
for i in range(1,9):
if m2 is None:
m2 = Basemap(llcrnrlat=min(lat2), llcrnrlon=min(lon2), urcrnrlat=max(lat2), urcrnrlon=max(lon2), resolution='f',
area_thresh=1)
ax.append(fig.add_subplot(3,3,i+1))
m2.drawcoastlines(linewidth=0.8)
#m[-1].shadedrelief()
im.append(m2.pcolormesh(lon2, lat2, um_data[i-1][0],vmin=0.0,vmax=5,cmap=colm))
tit.append('Ens. %s'%rom[i].upper())
ax[-1].set_title(tit[-1],fontsize=fontsize)
fig.subplots_adjust(right=0.99, bottom=0.15, top=0.95,left=0.01, hspace=0.2, wspace=0.0)
nplot+=1 ; _ = fig.savefig('Vid/Fig_%03i.png'%nplot, bbox_inches='tight', format='png', dpi=300)
#Plet the avg maps for comparison
fig = plt.figure(figsize=(18,12))
#sns.set_style("darkgrid", {'axes.grid':True, 'ticks':True})
UM_2 = UM133ens.groupby('t.day').sum(dim='t').mean(dim='day')
UM_1 = UM044ens.groupby('t.day').sum(dim='t').mean(dim='day')
obs = OBS.dataset[1].groupby('t.day').sum(dim='t').mean(dim='day')
colm = colmap2.Blues
colm.set_under('w', alpha=0)
tsteps = pd.DatetimeIndex(OBS.dataset[1].coords['t'].values).tz_localize(utc).tz_convert(timezone).to_pydatetime()
obs_data = np.ma.masked_less(obs.variables['lsrain'][:].values, -0.01)
data = [obs_data[1:-1]/6,
np.ma.masked_less(UM_1.variables['lsrain'][:].values, -0.01).mean(axis=0)[0]/6,
np.ma.masked_less(UM_2.variables['lsrain'][:].values, -0.01).mean(axis=0)[0]/6,
None,
np.ma.masked_less(UM_1.variables['lsrain'][:].values/6, -0.01).std(axis=0)[0],
np.ma.masked_less(UM_2.variables['lsrain'][:].values/6, -0.01).std(axis=0)[0]
]
data.append(None)
data.append(data[0] - data[1])
data.append(data[0] - data[2])
UM_1min = np.ma.masked_less(UM_1.variables['lsrain'][:].values, -0.01).min(axis=0)[0]
UM_2min = np.ma.masked_less(UM_2.variables['lsrain'][:].values, -0.01).min(axis=0)[0]
UM_1max = np.ma.masked_less(UM_1.variables['lsrain'][:].values, -0.01).max(axis=0)[0]
UM_2max = np.ma.masked_less(UM_2.variables['lsrain'][:].values, -0.01).max(axis=0)[0]
#data.append(None)
#data.append(UM_1max - UM_1min)
#data.append(UM_1max - UM_1min)
lat1, lon1 = obs.variables['latitude'][:,0], obs.variables['longitude'][0,:]
lat2, lon2 = UM_1.coords['lat'][:], UM_1.coords['lon'][:]
ax = [fig.add_subplot(3,3,1)]
#m[0].pcolormesh(topolon, topolat, ls.hillshade(topo[:], vert_exag=1), cmap='gray', ax=ax[0])
im = [m2.pcolormesh(lon1, lat1, obs_data/6,vmin=0.0,vmax=25/6,cmap=colm)]
m2.drawcoastlines()
#cbar_ax = fig.add_axes([0.14, 0.33, 0.74, 0.01])
#cbar=fig.colorbar(im[-1], cax=cbar_ax, orientation='horizontal')
#cbar.ax.tick_params(labelsize=24)
#cbar.set_label('Avg. Rain-rate [mm/d]',size=24)
tit = ['CPOL avg.', 'UM 0.44km ens avg.', 'UM 1.33km ens avg.',
None, 'UM 0.44 ens std.', 'UM 1.33km ens std.',
None, 'CPOL - UM 0.44 ens avg.' ,'CPOL - UM 1.33km ens avg.',
None, 'UM 0.44 ens range', 'UM 1.33km ens range']
colors = {1:colm, 2:colm, 4:colm, 5:colm,
7:colmap.GMT_polar_r, 8:colmap.GMT_polar_r,
10:colm, 11:colm}
Range = {1:(0.0, 25/6), 2:(0.0, 25/6), 4:(0,8/6), 5:(0,8/6), 7:(-20/6,20/6), 8:(-20/6,20/6), 10:(0,30/6), 11:(0,30/6)}
height = {2:0.682, 5:0.37, 8:0.044, 11:0.05}
ax[-1].set_title('CPOL avg.', fontsize=fontsize-2)
cbar_ax = [fig.add_axes([0.91, height[2], 0.01, 0.23])]
cbar = [fig.colorbar(im[-1], cax=cbar_ax[-1], orientation='vertical')]
cbar[-1].ax.tick_params(labelsize=fontsize-4)
cbar[-1].set_label('Rain-Rate [mm/d]',size=fontsize-4)
for i in range(1,len(data)):
if data[i] is not None:
ax.append(fig.add_subplot(3,3,i+1))
#m[-1].shadedrelief()
im.append(m2.pcolormesh(lon2, lat2, data[i],vmin=Range[i][0],vmax=Range[i][1],cmap=colors[i]))
m2.drawcoastlines()
ax[-1].set_title(tit[i],fontsize=fontsize-2)
if i in (5, 8, 11):
cbar_ax.append(fig.add_axes([0.91, height[i], 0.01, 0.23]))
cbar.append(fig.colorbar(im[-1], cax=cbar_ax[-1], orientation='vertical'))
cbar[-1].ax.tick_params(labelsize=fontsize-4)
cbar[-1].set_label('Rain-rate [mm/d]',size=fontsize-2)
fig.subplots_adjust(right=0.9, bottom=0.01, top=0.95,left=0.01, hspace=0.0, wspace=0)
nplot+=1 ; _ = fig.savefig('Vid/Fig_%i03.png'%nplot, bbox_inches='tight', format='png', dpi=300)
data, time, model = [], [], []
sns.set_style("darkgrid", {'axes.grid':True, 'ticks':True})
for i in range(len(ensembles)):
d1 = UM133.dataset[i].groupby('t.hour').mean(dim='t').mean(dim=('surface','lat','lon'))
d2 = UM044.dataset[i].groupby('t.hour').mean(dim='t').mean(dim=('surface','lat','lon'))
t1 = pd.Series(d1.variables['lsrain'].values, index=(d1.variables['hour'].values+9) % 24 ).sort_index(inplace=False)
t2 = pd.Series(d2.variables['lsrain'].values, index=(d2.variables['hour'].values+9) % 24 ).sort_index(inplace=False)
time += list(t1.index) + list(t2.index)
data += list(t1.values) + list(t2.values)
model += len(t1)*['UM 1.33km'] + len(t2)*['UM 0.44km']
cpol_1 = OBS.dataset[1].groupby('t.hour').mean(dim='t').mean(dim=('y','x'))
obs_S = pd.Series(cpol_1.variables['lsrain'].values, index=(cpol_1.variables['hour'].values+9) % 24 ).sort_index(inplace=False)
data += list(obs_S.values)
time += list(obs_S.index)
model += len(obs_S)*['Cpol']
plot_data = pd.DataFrame(dict(data=data, time=time, Type=model))
#nplot+=1 ; _ = fig.savefig('Vid/Fig_%03i.png'%nplot, bbox_inches='tight', format='png', dpi=300)
#mpld3.enable_notebook()
#Plot the maps
plt.close()
mpld3.disable_notebook()
sns.set_style("whitegrid", {'axes.grid':True, 'ticks':False})
um_1 = UM044ens.groupby('t.hour').mean(dim='t').mean(dim=('surface'))
um_2 = UM133ens.groupby('t.hour').mean(dim='t').mean(dim=('surface'))
obs = OBS.dataset[1].groupby('t.hour').mean(dim='t')
time = (um_1.hour + 9) % 24
data = [obs, um_2, um_1]
fig = plt.figure(figsize=(15,7))
lat1, lon1 = OBS.dataset[1].variables['latitude'].values[:,0], OBS.dataset[1].variables['longitude'].values[0,:]
lat2, lon2 = um_1.coords['lat'][:], um_1.coords['lon'][:]
fig.subplots_adjust(right=0.94, bottom=0.025, top=0.98,left=0.01, hspace=0, wspace=0)
first = True
val_max = 2
colm = colmap2.Blues
colm.set_under('w', alpha=0)
colm.set_bad('w', alpha=0)
#ax = [fig.add_subplot(2,3,1)]
mm, im = [], []
col = 2
for i in range(24):
#break
fname_1 = os.path.join(outdir, 'DiurnalCycle_%02i.png'%i)
if first:
for col in range(1, 7):
if col != 4:
ax1 = fig.add_subplot(2,3,col)
m2.drawcoastlines(linewidth=0.8)
if col == 1:
ax1.set_title('CPOL',fontsize=fontsize)
_ = im.append(m2.pcolormesh(lon1, lat1, obs['lsrain'].values[i], shading='gouraud', vmin=0.1, vmax=val_max, cmap=colm))
elif col == 2:
ax1.set_title('UM 1.33km mean', fontsize=fontsize)
_ = im.append(m2.pcolormesh(lon2, lat2, um_2['lsrain'].mean(dim='ens').values[i], shading='gouraud', vmin=0.1, vmax=val_max, cmap=colm))
elif col == 3:
ax1.set_title('UM 0.44km mean', fontsize=fontsize)
_ = im.append(m2.pcolormesh(lon2, lat2, um_1['lsrain'].mean(dim='ens').values[i], shading='gouraud', vmin=0.1, vmax=val_max, cmap=colm))
elif col == 5:
ax1.set_title('UM 1.33km std', fontsize=fontsize)
_ = im.append(m2.pcolormesh(lon2, lat2, um_2['lsrain'].std(dim='ens').values[i], shading='gouraud', vmin=0.1, vmax=val_max, cmap=colm))
elif col == 6:
ax1.set_title('UM 0.44km std', fontsize=fontsize)
_ = im.append(m2.pcolormesh(lon2, lat2, um_2['lsrain'].std(dim='ens').values[i], shading='gouraud', vmin=0.1, vmax=val_max, cmap=colm))
cbar_ax = fig.add_axes([0.14, 0.0, 0.74, 0.02])
cbar=fig.colorbar(im[-1], cax=cbar_ax, orientation='horizontal')
_ = cbar.ax.tick_params(labelsize=fontsize-4)
first = False
else:
###MEAN
im[0].set_array(obs['lsrain'].values[i].ravel())
im[1].set_array(um_2['lsrain'].mean(dim='ens').values[i].ravel())
im[2].set_array(um_1['lsrain'].mean(dim='ens').values[i].ravel())
###STD
im[-2].set_array(um_2['lsrain'].std(dim='ens').values[i].ravel())
im[-1].set_array(um_1['lsrain'].std(dim='ens').values[i].ravel())
cbar.set_label('Rain-Rate at %2i LT [mm/h]'%((i+9.5)%24),size=fontsize)
fig.savefig(fname_1, bbox_inches='tight', format='png', dpi=72)
#break
dest_dir = os.path.abspath('Vid')
make_mp4_from_frames(outdir, dest_dir, 'WeekOfHector-Diurnal-2.mp4', 4, glob='DiurnalCycle_??')
fig.clf()
plt.close()
day = (time[1:5] - 9) % 24
night = (time[5:5+4] - 9) % 24
obs_diff = obs['lsrain'][day].sum(dim='hour') - obs['lsrain'][night].sum('hour')
um_1_diff = um_1['lsrain'][day].sum(dim='hour') - um_1['lsrain'][night].sum('hour')
um_2_diff = um_2['lsrain'][day].sum(dim='hour') - um_2['lsrain'][night].sum('hour')
mask = np.ma.masked_invalid(obs['lsrain'].values) * 0 + 1
obs_max = (((np.argmax(obs['lsrain'].values * mask, axis=0)*mask[0] + 9.5 ) % 24)/2).round(0) * 2
um_1_max = ((um_1['lsrain'].argmax(dim='hour').mean(dim='ens').round(0) + 9.5) % 24 / 2).round(0) * 2
um_2_max = ((um_2['lsrain'].argmax(dim='hour').mean(dim='ens').round(0) + 9.5) % 24 / 2).round(0) * 2
colm = CubeYF_20.get_mpl_colormap(N=24, gamma=2.0)
#colm = qualitative.Paired[12].get_mpl_colormap(N=24, gamma=2.0)
colm.set_bad('w')
#m, ax = [], []
ax = []
data = [obs_max, um_1_max, um_2_max]
names = ['CPOL', 'UM 1.33km', 'UM 0.44km']
gs = gridspec.GridSpec(30, 3)
fig = plt.figure(figsize=(15,8))
fig.subplots_adjust(right=0.98, bottom=0.03, top=0.98,left=0.01, hspace=0, wspace=0)
cbar_ax = plt.subplot(gs[13, :])
for i in range(3):
ax.append(plt.subplot(gs[:13, i]))
if m2 is None:
m2 = Basemap(llcrnrlat=min(lat2), llcrnrlon=min(lon2), urcrnrlat=max(lat2), urcrnrlon=max(lon2), resolution='f',
area_thresh=1)
m2.drawcoastlines()
try:
im = m2.pcolormesh(lon1, lat1, data[i], vmin=0, vmax=24, cmap=colm)
except TypeError:
im = m2.pcolormesh(lon2, lat2, data[i], vmin=0, vmax=24, cmap=colm)
ax[-1].set_title(names[i], fontsize=fontsize)
cbar=fig.colorbar(im, cax=cbar_ax, orientation='horizontal',ticks=list(range(1,24,2)))
cbar.ax.tick_params(labelsize=fontsize-4)
#cbar.set_label('Local Time [h]', size=fontsize)
mpld3.disable_notebook()
# Plot avg diurnal cycle
#Get the minimum time period that is covered by all simulations and also observations
#fig = plt.figure(figsize=(15,8), dpi=72)
#ax = fig.add_subplot(111)
sns.set_style("darkgrid", {'axes.grid':True, 'ticks':True})
ax = plt.subplot(gs[18:, :])
####################
sns.lineplot(x=plot_data.time, y=plot_data.data, hue=plot_data.Type, ax=ax, lw=3, data=plot_data)
#ax.set_ylim(0.0,0.42)
ax.set_xlim(0,23)
ax.set_ylabel('Rain-Rate [mm/h]', fontsize=fontsize)
ax.set_xlabel('Local Time [h]', fontsize=fontsize)
ax.legend(loc=0, fontsize=fontsize-4)
ax.tick_params(labelsize=fontsize-4)
fig.subplots_adjust(bottom=0.145, top=0.99, right=0.99, left=0.1)
#ax.grid(color='r', linestyle='-', linewidth=0)
mpld3.disable_notebook()
nplot = 12; _ = fig.savefig('Vid/Fig_%03i.png'%nplot, bbox_inches='tight', format='png', dpi=300)
#Plot the pdfs
sns.set_style('darkgrid')
from mpl_toolkits.axes_grid.inset_locator import inset_axes
dist=lognorm([LogNorm.CPOL['shape']*3.4], loc = np.log(LogNorm.CPOL['scale']*1.9))#, shape # mu, sigma
med = np.median(um044vec.loc[um044vec>0].values)
X = np.linspace(0,65,200)
#b = 2.5
nbins = 100
def fitpdf(x, b):
a=med*100
return ((b/a) * (x/a)**(b-1)) / (1+(x/a)**b)**2
#popt, pcov = curve_fit(powerlaw, um044cnt['Rain-rate'].values[1:],
# um044cnt['count'].values[1:])
popt, pcov = PowerFit.CPOL['popt']/4, PowerFit.CPOL['pcov']
fig = plt.figure(figsize=(15,10))
ax = fig.add_subplot(111)
ax2 = ax.twinx()
#Define histogram data to plot
histdata = (um133vec.loc[um133vec>0], um044vec.loc[um044vec>0], obs_vec.loc[obs_vec>0])
#facecolors = ['fuchsia', 'firebrick', 'lightseagreen']
#n1_1, bins1_1, patches1_1 = ax2.hist(, bins=nbins, normed=1, facecolor='firebrick', alpha=0.25)
#n2, bins2, patches2 = ax2.hist(, bins=nbins, normed=1, facecolor='lightseagreen', alpha=0.25)
ax2.set_ylim(0,0.9)
ax2.set_yticks([])
#ax.plot(np.linspace(0.01, 100,200), fitpdf(np.linspace(0.01,100,200),popt[0]), lw=4, color='navy',
# label='Log-logistic $\\beta = %02.2e$'%popt[0])
#ax.plot(cnts.UM044['Rain-rate'].values, cnts.UM044['counts'].values,'-', label='UM 0.44 km ens',lw=4, color='fuchsia')
#ax.plot(cnts.UM133['Rain-rate'].values, cnts.UM133['counts'].values,'--', label='UM 1.33 km ens',lw=4, color='fuchsia')
#ax.plot(cnts.CPOL['Rain-rate'].values, cnts.CPOL['counts'].values,'-', label='CPOL',lw=4, color='lightseagreen')
#ax.scatter(um044cnt['Rain-rate'].values, um044cnt['count'].values, c='red', s=125, alpha=0.24)
#ax.scatter(obscnt['Rain-rate'].values, obscnt['count'].values, c='purple', s=125, alpha=0.24)
#ax.plot(np.linspace(0, 50,200), um044pdf, lw=4, label='UM 0.44km ens', color='lightseagreen')
#ax.plot(np.linspace(0, 50,200), obspdf, lw=4, label='CPOL', color='fuchsia')
#ax.plot(np.linspace(0.0,50,200), dist.pdf(np.linspace(0.0, 50,200)), color='firebrick', alpha=0.8,
# lw=4, label='Log-norm ($\\mu=%02.2f, \\sigma=%02.2f$)'%(scale, shape))
n1, bins1, patches1 = ax.hist(histdata, bins=nbins, normed=1, alpha=0.55)
#fitted_pl.plot_pdf(ax=ax,color='k',lw=4)
#ax.plot(np.linspace(0.1, 100), um044dens_aavg(np.linspace(0.1,100)), lw=4, label='UM 0.44km area-avg')
#ax.plot(np.linspace(0.1, 100), um133dens_aavg(np.linspace(0.1,100)), lw=4, label='UM 1.33km area-avg')
ax.tick_params(labelsize=fontsize-4)
#ax.set_ylim(0.,.25)
ax.set_xlim(0.,20)
ax.set_ylabel('Density', fontsize=fontsize)
ax.set_xlabel('Rain-Rate [mm/h]', fontsize=fontsize)
axin = inset_axes(ax, width = "80%", height = "80%", loc=1)
axin.plot(cnts.UM133['Rain-rate'].values, cnts.UM133['counts'].values,'-', label='UM 1.33 km ens',lw=2)
axin.plot(cnts.UM044['Rain-rate'].values, cnts.UM044['counts'].values,'-', label='UM 0.44 km ens',lw=2)
axin.plot(cnts.CPOL['Rain-rate'].values, cnts.CPOL['counts'].values,'-', label='CPOL',lw=2)
#axin.plot(np.linspace(0.1,80,200), dist.pdf(np.linspace(0.1, 80,200)), lw=1, alpha=0.8,
# label='Log-norm ($\\mu=%02.2f, \\sigma=%02.2f$)'%(scale, shape))
#axin.set_xlabel('$\\log$(Rain-rate)', fontsize=fontsize-2)
#axin.plot(X, um044pdf,'-', label='UM 0.44 km ens',lw=4, color='fuchsia')
#axin.plot(X, obspdf,'-', label='CPOL',lw=4, color='lightseagreen')
axin.set_yscale('log')
axin.set_xscale('log')
#axin.set_xticks([])
#axin.set_yticks([])
axin.tick_params(labelsize=fontsize-4)
axin.set_xlim(0.1, 100)
axin.set_ylim(0.000023, 1)
ann = axin.annotate("(b)", (0.95,0.93), xycoords="axes fraction", fontsize=fontsize)
axin.legend(loc=3, fontsize=fontsize-4)
#n.set_ylabel('$\\log$(Density)', fontsize=fontsize-2)
fig.subplots_adjust(bottom=0.13, right=0.98, top=0.97)
nplot+=1 ; _ = fig.savefig('Vid/Fig_%03i.png'%nplot, bbox_inches='tight', format='png', dpi=300)
#mpld3.enable_notebook()
Now get the percentiles for the ensemble simulation and compare it against the observations
#Plot the percentages
fig = plt.figure(figsize=(15,10))
ax = fig.add_subplot(111)
plt.subplots_adjust(bottom=0.05, right=0.9, top=0.4)
ax.plot(PERC.index, PERC['UM 1.33km'].values, label='UM 1.33km',lw=3)
ax.plot(PERC.index, PERC['UM 0.44km'].values, label='UM 0.44km', lw=3)
ax.plot(PERC.index, PERC['Obs'].values, label='CPOL',lw=2)
ax.set_xlim(1,100)
#ax.set_ylim(10,100)
#ax.set_yscale('log')
#ax.set_xscale('log')
ax.set_xlabel('Percentile []', fontsize=fontsize)
ax.set_ylabel('Rain-rate [mm/h]', fontsize=fontsize)
ax.tick_params(labelsize=fontsize-4)
axin = inset_axes(ax, width = "80%", height = "60%", loc=9)
axin.plot(PERC.index[1:], PERC['UM 1.33km'].values[1:], label='UM 1.33km',lw=3)
axin.plot(PERC.index[1:], PERC['UM 0.44km'].values[1:], label='UM 0.44km', lw=3)
axin.plot(PERC.index[1:], PERC['Obs'].values[1:], label='CPOL',lw=2)
#axin.set_xticks([])
#axin.set_yticks([])
axin.tick_params(labelsize=fontsize-4)
axin.set_xlim(50, 100)
axin.set_ylim(.25,100)
#axin.set_xticks([20, 200, 500])
axin.legend(loc=0, fontsize=fontsize-4)
#axin.get_xaxis().set_major_formatter(mpl.ticker.ScalarFormatter(useMathText=False))
#axin.get_xaxis().get_major_formatter().labelOnlyBase = True
#axin.set_xticklabels([], rotation=0)
ann = axin.annotate("(b)", (0.95,0.93), xycoords="axes fraction", fontsize=fontsize)
#axin.set_xticks(range(50,100,10), minor=False)
axin.set_yscale('log')
axin.set_xscale('log')
fig.subplots_adjust(bottom=0.13, right=0.98, top=0.97)
nplot+=1 ; _ = fig.savefig('Vid/Fig_%03i.png'%nplot, bbox_inches='tight', format='png', dpi=300)
nplot
15
mpl.ticker.ScalarFormatter?